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Abstract: The development of a pressure-dependent constitutive model with combined multilinear 
kinematic and isotropic hardening is presented. The constitutive model is developed using the 
ABAQUS user material subroutine (UMA T). First the pressure-dependent plasticity model is 
derived. Following this, the combined bilinear and combined multilinear hardening equations are 
developed for von Mises plasticity theory. The hardening rule equations are then modified to 
include pressure dependency. The method for implementing the new constitutive model into 
ABAQUS is given. Finally, verification of the UMAT is presented. 
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1 . Introduction 


Classical metal plasticity theory assumes that hydrostatic stress has a negligible effect on the yield 
behavior of metals (Bridgman, 1947; Hill, 1950). Spitzig (1975, 1976,1984) and Richmond 
(1980) conducted experiments demonstrating that yielding in metals could be approximated as a 
linear function of hydrostatic stress. Recent reexaminations of classical theory by Wilson (2002) 
and Allen (2000,2002) have revealed a significant effect of hydrostatic stress on the yield behavior 
of various metals. An ABAQUS user material subroutine (UMAT) was developed to investigate 
the effect of hydrostatic stress on low cycle fatigue of metals. The UMAT is based on the 
Drucker-Prager (1952) yield theory and incorporates combined multilinear kinematic and isotropic 
hardening. A summary of the development of the UMAT is presented here. See Allen (2002) for 
more details on the UMAT development including the corresponding UMAT FORTRAN code. 
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2. Pressure-dependent plasticity model 


This section details the development of the governing equations for a pressure-dependent plasticity 
model including numerical methods for the integration of the constitutive equations. The 
equations are derived using the method developed by Aravas (1987,1996). The tensor 
components are given with respect to a cartesian coordinate system with indices ranging from 1 to 
3. 

2.1 Yield function 

The chosen pressure-dependent yield function is the Drucker-Prager yield function. The Drucker- 
Prager yield function is written in ABAQUS (ABAQUS, 1995) notation as 

f =t- ptan0-d = 0, (1) 

where t is a pseudo-effective stress , 6 is the slope of the linear yield surface in the p-t stress plane, 
p is the hydrostatic pressure, d is the effective cohesion of the material. In metal plasticity, d is 
equivalent to the current yield stress, <t vs . In general, er„ is a function of the equivalent plastic 
strain, , which is given by the equation 


£ p‘ = j- £ P‘ £ P l 

eq V 3 ‘J ’ 

where £y is the plastic strain ( £? 1 = £y -£y ). 

The variables used by the linear Drucker-Prager yield function are shown graphically in Figure 1. 
The flow potential, g, for the linear Drucker-Prager model is defined as 

g = t - ptanyr , 


( 2 ) 


(3) 


where y/is the dilation angle in the p-t plane. Setting y/= £? results in associated flow. Therefore, 
the original Drucker-Prager model is available by setting y/= 9 and r = 1, which leads to 

t = V 3 -/ 2 ■ Substituting p = -^J l into equation (1) gives 

/ = V^2 +-^I\tand -d , 

where /, is the first stress invariant. To conveniently compare ABAQUS Drucker-Prager 
material property variables with those used by Richmond, et al. (1980) in their material testing let 

a = — tan# . 

3 

Finally, Equation (4) can be written as 

/ 


(4) 


( 5 ) 
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f = °eff - 3fl/> - <T VS (eQ ) = 0 , 

where cr^ is the effective stress for von Mises plasticity theory. 


( 6 ) 


2.2 Flow Rule 

The rate form of the associated flow rule is written as 


d(Ty 


The flow rule is written in a more general form as 


• pi • — 1 • T 

4 = + jV’ 


where 


ft:: = 


1 S: 

IJ 


ij 2 <v’ 


4 =# 




da, 


£ P =~9 


eff 

& 

dp ’ 


where Sy is the deviatoric stress tensor and £ q and £ p are the distortional and volumetric portions 

of the plastic strain rate, respectively. Eliminating (j) from Equations (10) and (11) gives 

. df . df . 

r» — + — = 0, 


q dp p da, 


eff 


and Equation (12) can be rewritten as 


4 (-3 a)+d p =0. 


(7) 

( 8 ) 

(9) 

( 10 ) 

(ID 


( 12 ) 


(13) 


2.3 State variable equations 

For the case of the linear Drucker-Prager model, only one state variable, f 1 is required and is 
defined as 

r 1 s e * . (14) 

The plastic work equation is used to derive a usable form of the state variable equation, which can 
be written as 
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( 15 ) 


£ Pl =■ 
b eq 


- P £p+(T eff £ q 


■>>sV / 


2.4 Numerical integration 

The backward (implicit) Euler method was chosen for the integration of the elastoplastic 
equations. Ortiz and Popov (1985) found that the backward Euler method is very accurate for 
strain increments several times the size of the yield surface in strain space and that the method is 
unconditionally stable. 

Assuming isotropic linear elasticity, one can write 

<*ij = 1 V£y +**«!*. ( 16 ) 


where a t] is the stress tensor, p is the shear modulus, K is the bulk modulus, and I ; . is the second- 

order identity tensor. Using the backward Euler method to integrate equations (16), (8), (13) and 
(15) leads to the following incremental forms 


a u L =°f - KA£ P hj’ 

(17) 

\£y - A£ q Tly | B+J + ~ teplg , 

(18) 

A£ q {-3a)+A£ p =0, 

(19) 

a ys / 

(20) 


The “n+1” subscript denotes values at the end of the increment at time, t| n+] =t\ n + At , and the 
“pr” superscript indicates a predicted value based on the purely elastic solution. 

Summarizing the primary equations and dropping the “n+1” subscript for brevity, one can write 


a^-3ap~a ys (£^)=0. 

(21) 

A£ q (-3a)+A£ p =0, 

(22) 


(23) 

p = p pr +KA£ p , 

(24) 


and 
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( 25 ) 


* ' 

This forms a set of five nonlinear equations that is solved for p, ASp, Ae q , and A . Solving 

Equations (21) through (25) for the five unknowns gives 

p pr p + Kaa£ - Kaajtfj , ) 


P = 


a eff =' 




3Ka 2 + p 

3 ap pr p + 3 Ka^a 2 + per ys [e^) 
3Ka 2 + p 


a|rr 


A £„ = 


3Ka 2 + p 

~pr — l~pt\ 

“eff -“ys\*e q r 2a P 

9Ka 2 + 3p 


and 


Ae pl = 

eg 


° P eff - a ys( £ ^)~ 3a P Pr 


9 Ka 2 + 3p 


In general, Newton iteration is used to solve Equation (30) for A and a ys [e p q ), and then 

Equations (26) through (29) are solved by direct substitution. Strain increments and stresses are 
then determined by 

1 


A£ y =-A^iy +h£ g h ij 


and 


a ij Plij + j a eff n ij 

thereby completing the solution of the elastoplastic equations. 

In an implicit finite element code, such as ABAQUS, the equilibrium equations are written at the 
end of the increment resulting in a set of nonlinear equations in terms of the nodal unknowns. If a 
Newton scheme is used to solve the global nonlinear equations, the Jacobian (tangent stiffness) 
must be calculated (Aravas, 1996). The Jacobian, J is defined as 


J = 


da\ 

H 


(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 


(33) 


2004 ABAQUS Users ' Conference 


5 


The accuracy of Jacobian effects convergence, and, therefore, errors in the Jacobian formulation 
may result in analyses that require more iterations or, in some cases, diverge (ABAQUS, 2001). 
The Jacobian for von Mises plasticity was employed in this research. 


3. Hardening models 


Next, the development of the combined kinematic and isotropic hardening models is presented. 
For simplicity, die hardening equations for the classical von Mises theory along with a bilinear 
hardening theory are developed first. The derivation of the bilinear hardening equations closely 
parallels die work of Taylor and Flanagan (1998). These equations are then modified to form the 
more complex multilinear hardening equations. Finally, the bilinear and multilinear equations are 
modified to include the Drucker-Prager pressure dependency. To begin, several terms that are 
common to each of the hardening models are defined. 

3.1 Basic definitions 


The center of the yield surface in deviatoric stress space (the backstress) is defined by the tensor 
Oij. The deviatoric stress tensor, S tJ is defined as 

s ij - a ij ~~ l i ■ 

The stress difference, 4 is defined by 

£ij ~ '-* !/ ~ ®| i ■ 

The magnitude of the stress difference, R, is then written as 


(34) 

(35) 

(36) 


The geometric relationship between gy, ay, and Sy leads to taking the magnitude of the stress 
difference as a radius, hence the symbol R. The von Mises yield surface in terms of the stress 
difference is 


f = ±4 iJ S j =k 2 , (37) 

where it is the yield strength in pure shear, and the effective stress is 

< 38 > 

Combining Equations (36) and (38), one defines 

R 


R as 
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( 40 ) 


The normal to the yield surface, Q tJ is then determined from Equation (37) as 

’ l^/^sl R 

Assuming a normality condition, the plastic part of the strain is 

£$ l =rQij, (41) 

where y is a scalar multiplier that must be determined (Taylor, 1998). 

3.2 Combined bilinear hardening 


Assuming a linear combination of the two hardening types, a scalar parameter, p, can be defined 
which determines the amount of each type of hardening with 0 < /? < 1 . A value of P= 1 
indicates only isotropic hardening, and a value of P = 0 indicates only kinematic hardening. 


Equations for R and ay are written as 



and 

d ij = ^HrQ ij (i-p), 


(42) 


(43) 


where H is the slope of the equivalent stress versus the equivalent plastic strain as shown in Figure 
2. The consistency equation is formed as 


Qy4y-R 


or 



Using the additive strain rate decomposition and the elastic stress rate with Equation (45) and 
taking the tensor product with the normal, £?/, gives 


Qij^ij -iQijDykiQki -Qy 


-^Hrii-p) 


Qn 



(44) 


(45) 


(46) 


Solving for /gives 
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1 + 


H_ 

in 


Qij £ ij ■ 


The incremental forms of the governing equations are 

”• Lr<L- Ar 2 f,e 'J- 




and 


a ‘j\ n+ r a ul + ^ HA rQ 0 (i-/J) 

where, as before, the subscripts “n” and “n+1” refer to the beginning and end of a time step, 
respectively. 

An incremental form of the consistency condition is also required and is written as 

CC;: I +R I O.. = S.-J 

»l«+i '"+i v y| n+1 

Substituting Equations (48) through (50) into the consistency condition of (51) gives 


a U L +\ha y p 


Qij = S r Ay 2 MQiJ . 


3 "j L 3 

Taking the tensor product of both sides of the previous equation with Q tj and solving for Ay gives 


Ay=- 


1 




i + " 
3/T 


n 


p r \ 


«n+l 


R l • 


(47) 


(48) 

(49) 

(50) 


(51) 


(52) 


(53) 


Equation (53) demonstrates that the plastic strain increment is proportional to the magnitude of the 
distance of the elastic predictor stress past the yield surface Having now determined Ay from 
Equation (53) Equations (48) through (50) can be solved. In addition, one also computes 


Ae[i = QijAy 


and 


A< 


‘ill ir ’ 


which completes the bilinear combined hardening incremental solution (Taylor, 1998). 


(54) 


( 55 ) 
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3.3 


Combined multilinear hardening 


The equations for the multilinear hardening case are similar to those for the bilinear case, but an 
added level of complexity is required because the slope of the equivalent stress versus equivalent 
plastic strain, H is no longer constant. This complication is clearly seen by examining Figure 2, 
where for a given plastic strain increment, //may take on several values, //,. 

Recall that for bilinear hardening Aa^ was written as 


Aa^-HAs* 


Similarly, for the multilinear hardening case Acr^ can be written as 

A(J eff = h( e £ Y £ £ = *(< )Jj A r . 

where h(s£j ) is a hardening function. Substituting k{e ^ ) for H in Equations 42 and 43 and 


integrating gives 


and 




4 a > [« MO] G «L ■ 

Rewriting Equation (5 1) in terms of Equations (58) and (59) (where /?| n+J = /?| n + AR , etc.) gives 

- a 'i ={ +s i. + iIK £ «l 

Multiplying both sides of the previous equation by Qy | ^ leads to 


& 


+t 


4 <MO 




( 56 ) 


(57) 


(58) 


(59) 


(60) 


(61) 


which is a nonlinear equation that can be solved for Ae ^ . Once Ae^ is known, the remaining 

incremental elastoplastic equations can be solved following the method demonstrated for the 
bilinear hardening case. 
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3.4 Pressure-dependent combined hardening 

To conveniently compare the pressure-independent combined hardening equations with the 
previously derived Drucker-Prager elastoplastic equations, the following substitutions are made: 


and for the multilinear case 


and 




Making the previous substitutions into Equation (53) gives 

a *- a A, 


as p ‘ = 

** 3 fj + H 


Recall the equation for pressure-dependent equivalent plastic strain. Equation (30) 

29 9Ka 2 +3fj 

Substituting a = 0 for pressure-independent plasticity and 

^k'L. =<r #\ m +HAe $ 

for bilinear hardening results in repeating equation (67). 

In a similar manner for the multilinear hardening case, substituting Equations (62) through (66) 
into Equation (61) and solving for gives 




Ge ff^\„ + \ +<Te ff\n a A n 


( 62 ) 

(63) 

(64) 

(65) 

( 66 ) 

(67) 


( 68 ) 


(69) 


l/H-l 

3 // 


(70) 
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For the case of pure isotropic hardening 



and Equation (70) can be written as 


A<=- 


'S-^k L 

3 // 


which is identical to Equation (68) for pressure-independent (a = 0) plasticity. Therefore, adding 
the terms for a ^ | and a vs |^ to the numerator of Equation (68) results in the proper form of the 

equation for Ae^ for pressure-dependent combined multilinear hardening, which is written in its 
final form as 


= 


T Pr 

eff 


~ a ys(^\ n 


+ <y. 


eff 


- a 




3 ap 


Pr 


9 An 2 + 3p 


( 71 ) 


(72) 


4. Constitutive model programming 


AJBAQUS provides a useful user subroutine interface called UMAT that allows one to define 
complex or novel constitutive models that are not available with the built-in ABAQUS material 
models. UMAT’s are written in FORTRAN code, and these FORTRAN subroutines are linked, 
compiled, and used by ABAQUS in the finite element analysis. 

Two UMAT’s were developed for this research. The first UMAT embodied a Drucker-Prager 
constitutive model with combined bilinear hardening. This subroutine was developed as a test 
case to confirm the validity of the pressure-dependent and combined hardening equations. The 
second UMAT embodied a Drucker-Prager constitutive model with combined multilinear 
hardening. This subroutine was developed for use in the low cycle fatigue metal plasticity 
analyses. Both subroutines were written in FORTRAN and are included in Allen (2002). 

5. UMAT program verification 


An axisymmetric finite element model of a notched round bar (NRB) geometry with the notch root 
radius, p, equal to 0.040 in. was created to simulate NRB tensile and low cycle fatigue tests. The 
NRB geometry was chosen to represent a specimen with a high hydrostatic stress influence (Allen, 
2000). The NRB finite element model was composed of approximately 500 Q4 axisymmetric 
elements (type CAX4 in ABAQUS). Two planes of symmetry were utilized as illustrated in 
Figure 3. A uniform displacement was applied to the top nodes of the FEM’s for the loading 
boundary condition 

The user defined constitutive model (UMAT) was compared to several of the ABAQUS built-in 
material models to test the accuracy of the elastoplastic equations and the corresponding 
FORTRAN code. Aluminum 2024-T851 material properties reported by Allen (2002) were used 
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in the NRB finite element model. The Drucker-Prager constant, a was 0 for von Mises plasticity 
and by trial and error was assumed to be 0.041 for pressure-dependent plasticity. 

The first test case was the monotonic loading of the NRB FEM using the multilinear isotropic von 
Mises and multilinear isotropic Drucker-Prager built-in constitutive models in ABAQUS. For this 
case, p = 1 (pure isotropic hardening) for the combined multilinear hardening UMAT. The results 
of this test are shown in Figure 4, and clearly show that the global responses of the UMAT and the 
built-in ABAQUS models are identical. 

The next test case was the cyclic loading of the NRB FEM using the bilinear kinematic von Mises 
model in ABAQUS. ABAQUS does not have a built-in model for the Drucker-Prager model with 
kinematic hardening or for the von Mises model with multilinear kinematic hardening. For this 
case, only the first and last values in the cr- e pl data table were used with the combined multilinear 
hardening UMAT to simulate bilinear hardening, and p = 0 (pure kinematic hardening). The 
results of this test are shown in Figure 5, which once again show that the global responses of the 
UMAT and the built-in ABAQUS model are identical. The Drucker-Prager UMAT curve is 
included in Figure 5 for comparison and follows the expected trend. 

The final test case was to examine the response of the combined multilinear hardening UMAT for 
varying p values. Cyclic finite element analyses with the NRB were conducted with /? = 0, 0.5, 
and 1. The results from the first cycle are shown in Figure 6. The three FEM’s have practically 
identical load-gage displacement responses until the negative loads are reached on the first 
unloading cycle. From this point on, the three solutions diverge in the expected manner with the p 
= 0.5 solution falling between the pure isotropic and pure kinematic curves. 
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7. Figures 



Figure 1. Linear Drucker-Prager model: yield surface and flow direction in the p-t 
plane (Adapted from (ABAQUS, 1995)) 
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Figure 2. Illustration of the relationship between yield stress and equivalent plastic 
strain for the (a) bilinear hardening and (b) multilinear hardening case 



Figure 3. Schematic of axisymmetric model of a notched round bar specimen with 
p= 0.040 in. utilizing two planes of symmetry 
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Figure 4. Comparison of the built-in ABAQUS models with multilinear isotropic 
hardening and the combined multilinear hardening UMAT 
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Figure 5. Comparison of the built-in ABAQUS von Mises model with bilinear 
kinematic hardening and the combined multilinear hardening UMAT 
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Figure 6. Comparison of Drucker-Prager combined multilinear hardening UMAT 

solutions with different ^Values 

8. Acknowledgement 


G. 0 


k . + *b 
** o 
c o c x x 


+ ;/ 


/ 


0.5 in. Gage Length 

Drucker-Prager UMAT. 0 = 1 .0 
Drucker-Prager UMAT, 0 = 0.5 
Drucker-Prager UMAT, 0 = 0.0 


We would like to thank several people at Marshall Space Flight Center for their assistance, 
guidance, and advice. These people include Dr. Greg Swanson, Jeff Rayburn, Dr. Preston McGill, 
and Doug Wells. In addition, we would like to express thanks to Bill Scherzinger and Kevin 
Brown at Sandia National Laboratories for their technical assistance. Funding for this research 
was provided by the NASA Graduate Student Researchers Program (GSRP). Support was also 
provided by the Center for Manufacturing Research and the Mechanical Engineering Department 
at Tennessee Technological University. 


16 


2004 ABAQUS Users ’ Conference 



